function y = func_Wmn(omega,chi,m,n,CONSTS)

    c     = CONSTS.c;
    a     = CONSTS.a;
    d = CONSTS.d;
    eps_a = CONSTS.eps_a;
    wH    = CONSTS.omega_H;
    wp    = CONSTS.omega_p;
    
    k0 = omega./c;
    Z0 = 4*pi/c;
    epsVSw = (1 + (wp^2)./(wH^2 - omega.^2));
    etaVSw = (1 - (wp^2)./(omega.^2));
    ksi_square = (1/4)*(abs(epsVSw.*etaVSw)+eps_a^2);
    hi0 = atan((eps_a)./(sqrt(abs(epsVSw.*etaVSw))));
    h_mnVSw = (1./(a)).*sqrt(abs(epsVSw./etaVSw)).*(hi0 + (pi/2)*(2*n+abs(m)+1/2).*(sign(epsVSw)));
    
    y = (1/(8*pi))*Z0*k0*(chi.^2).*(epsVSw./h_mnVSw).*((m^2)./((m^4)+(a^4*k0.^4.*ksi_square)-(a^2*k0.^2.*m^2.*eps_a))).*...
        ((sin(h_mnVSw.*d))./(h_mnVSw.*d)).*(besselj(0,(h_mnVSw.*d)));    
end

